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ABSTRACT 

Using  analytical  and  numerical  techniques  we  compare  the  effect  of  three 
different  time-differencing  methods  on  the  propagation  of  shallow-water  gravity 
waves.  We  compare  the  numerical  dispersion  of  both  phase  and  group  velocities 
caused  by  a  centered-in-time  explicit  formulation,  a  centered-in-time  implicit 
method,  and  a  trapezoidal -in- time  implicit  treatment.  All  three  techniques  slow 
the  waves  relative  to  their  theoretical  phase  and  group  velocities.  For  each 
method,  the  slowing  Increases  with  decreasing  spatial  resolution  and  with  decreasing 
temporal  resolution.  For  At  s/giT/d  <  1,  the  slowing  increases  as  we  switch  from 
explicit,  to  trapezoidal  implicit,  to  centered  implicit  methods.  In  this  formula, 
x/gH1  =  theoretical  shallow  water  gravity  wave  speed.  At  =  time  increment,  and 
d  =  smallest  space  increment  between  like  variable  grid  points.  For  At,/gH/d  >  1, 
the  trapezoidal  implicit  scheme  again  outperforms  the  centered  implicit  scheme.. 
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I .  Introduction 

Implicit  time  differencing  in  a  finite  difference  approach  to  an  initial 
value  problem  indicates  that  the  solution  for  a  given  variable  at  a  given  time 
level  requires  knowledge  of  variables  at  the  same  time  level  in  the  forcing 
functions.  In  an  explicit  technique,  only  knowledge  of  variables  at  previous 
time  levels  is  required.  This  study  discusses  a  semi-implicit  scheme  for  the 
equations  of  motion  that  treats  gravity  waves  implicitly,  but  terms  for  other 
physical  processes  explicitly.  This  technique  numerically  retards  the  phase 
speed  of  the  shallow-water  gravity  waves.  This  slowing  then  allows  the  use  of 
large  time-steps,  which  would  violate  the  Courant  limitation  affecting  explicit 
schemes  (Mesinger  and  Arakawa,  1976). 

Kwizak  and  Robert  (1971)  pioneered  the  use  of  the  semi-implicit  technique 
in  an  atmospheric  grid  point  model.  O'Brien  and  Hurlburt  (1972)  and  Hamilton 
and  Rattray  (1978),  among  others,  applied  this  method  in  coastal  upwelling 
problems.  Hamilton  (1976)  and  Wang  and  Kravitz  (1980)  used  this  form  of  time 
differencing  in  applications  to  estuarine  circulation.  All  these  applications 
include  possible  effects  of  gravity  waves,  but  only  Wang  and  Kravitz  (1980) 
mentioned  any  effect  of  time-differencing  Induced  numerical  dispersion  on  their 
results.  They  note  the  retarding  effect  of  the  semi-implicit  scheme  by  stating 
that  it  is  "apparently  the  result  of  damping  of  surface  gravity  waves  by  the 
semi-implicit  method."  They  further  discuss  the  apparent  limitations  of  the 
method  In  transient  conditions  when  seiching  becomes  considerable.  One  can 
generalize  the  limitations  of  the  semi-implicit  approach  to  any  situation  where 
shallow  water  gravity  waves  are  significant.  Including  ocean  forecasting  in 
coastal  areas  where  transient  forcing  due  to  winds  or  tsunamis  may  arise. 

Following  the  methods  of  Mesinger  and  Arakawa  (1976)  and  Grotjahn  and  O'Brien 
(1976),  we  demonstrate  the  actual  effect  of  the  semi-implicit  technique  on  the 
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phase  and  group  velocities  of  shallow-water  gravity  waves.  We  compare  the 
analytic  results  with  a  centered  space  and  time  explicit  technique,  the  centered 
space  and  time  implicit  method  of  Kwizak  and  Robert  (1971)  and  the  centered  In 
space  trapezoidal  in  time  method  used  by  Wang  and  Kravltz  (1980).  We  Include 
the  additional  slowing  effects  of  the  space  differencing  to  make  our  results  as 
comparable  as  possible  to  numerical  simulation  applications.  For  clarity  In  the 
analysis,  we  shall  ignore  all  terms  but  those  generating  shallow-water  gravity 
waves . 

2.  Linearized  Shallow  Water  Gravity  Wave  Equations 

The  one-dimensional,  linearized,  shallow-water  gravity  wave  equations  may 
be  written 

ut  =  (2.1) 

nt  -  -Hux  (2.2) 

where  u  =  horizontal  velocity,  g  *  gravitational  acceleration,  n  =  surface  height 
deviation,  and  H  *  total  depth;  subscripts  indicate  differentiation  by  the  sub¬ 
scripted  variable. 

We  assume  wave  form  solutions 


u(x,t)  =  RelGe1^-"4'] 

(2.3) 

(2.4) 

with  k  =  radian  wave  number  and  w  «  radian  wave  frequency.  We  then  substitute 
(2.3)  and  (2.4)  into  (2.1)  and  (2.2)  treating  H  as  a  constant.  These  solutions 
satisfy  (2.1)  and  (2.2)  for  the  frequency 

u)=  +  k/gH.  (2.5) 

We  calculate  phase  velocity,  C,  from 

C  *  *  +  /giT  (2.6) 
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and  group  velocity,  Cg,  from 


dw  _  .  yi^? 
HIT  ‘  -  gH  * 


These  results  tell  us  that  the  solutions  to  (2.1)  and  (2.2)  are  nondlspersive 
waves  with  equal  phase  and  group  speeds  that  depend  only  on  gravitational 
acceleration  and  water  depth.  Recall  that  the  group  speed  represents  the  energy 
propagation  speed. 

We  will  demonstrate  that  solutions  to  (2.1)  and  (2.2)  by  finite  difference 
techniques  introduce  numerical  dispersion  into  the  solutions  and  that  this 
dispersion  is  highly  technique  dependent. 

3.  Explicit  Treatment.  Centered  in  Space  and  Time 


We  write  the  explicit,  centered  in  space  and  time,  finite  difference  form  of 
(2.1)  and  (2.2)  as 


(n+l)At  (n-l)At 
ujd*  ujd* 


i+l)d*  - 


(n+l)At  (n-l)At 
njd*  nid* 


nAt  nAt 


(3.1) 


(3.2) 


where  At  *  finite  time  increment,  d*  =  finite  space  increment  between  grid  points*, 
superscripts  represent  position  in  time,  and  subscripts  indicate  position  in  space. 


We  now  write  solutions  (2.3)  and  (2.4)  as 

u(jd*,nAt)  =  Relue1^*1*  "  unAt)] 
n(jd*,nAt)  *  Re(ne1{kJd*  '  unAt)]. 

Substituting  (3.3),  (3.4)  into  (3.1),  (3.2)  yields  solutions  if 
w*  *  +  sin"*  [^>^7  sin(kd*)] 


(3.3) 

(3.4) 


(3.5) 


*Note  d*  *  %d  where  d  Is  distance  between  like  labeled  grid  points  on  Arakawa  grid 
c.  (Mesinger  and  Arakawa,  1976). 
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This  result  then  implies  a  dispersive  wave  with  phase  speed,  C*,  and  group  speed. 


C*  =  =  -  kit  sin"lf3^  s-in(kd*) 3  (3.6) 

C|  =  ^  =  +  ^H*  cos(kd*)  [1  -  sin^kd*)]’5*  .  (3.7) 

We  now  introduce  the  effects  of  a  time-centered  implicit  formulation. 


4.  Implicit  Treatment,  Centered  in  Space  and  Time 

Kwizak  and  Robert  (1971)  use  the  implicit,  centered  in  space  and  time, 
finite  difference  form  of  (2.1)  and  (2.2),  which  we  write. 


,(n+l)At  ..(n-l)At 


(n-l)At  (n-l)At  (n+l)At  (n+l)At 


jd*  ~  ujd*  _  22.  rn(-i+l)d*  ~  n(j-l)d*1  n(,j+l)d*  ~  n(,i-l)d*  , 

2At  ”  2  1  fd*  2d*  1 


(4.1) 


(n+l)At  (n-l)At 
njd*  “  njd* 

2At 


-H 

~r 


(n-l)At  (n-l)At 
jU(j+l)d*  u(j-l)d* 


2d* 


u 


(n+l)At  „(n+l)At 
(j*l)d*;U(j-l)d* 


2d* 


■] 


(4.2) 


Substituting  the  finite  difference  solutions  (3.3)  and  (3.4)  into  (4.1)  and  (4.2) 
now  yields  a  different  numerically  generated  dispersion.  The  radian  frequency 
now  involves  the  inverse  tangent; 

w*  =  -  "St  tan_1  sin(kd*)l .  (4.3) 

The  phase  and  group  speeds  tlien  are 

C*  =  +  tan  ^  [^4  sin(kd*))  (4.4) 


and 

CJ  -  =  i^iH  cos(kd*)[l  +  ($r  v®,)2sin2(kd*)]‘1.  (4.5) 

Note  that  no  information  at  the  nth  time  level  appears  above.  The  significance  of 
this  becomes  apparent  as  we  look  at  a  trapezoidal  Implicit  scheme  similar  to  that 
used  by  Wang  and  Kravitz  (1980). 
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5.  Implicit  Treatment,  Centered  in  Space,  Trapezoidal  in  Time 

We  write  the  implicit,  centered  In  space,  trapezoidal  in  time,  finite 
difference  form  of  (2.1)  and  (2.2)  as 


(n+l)At  „nAt 
Jjd* 


At 


nAt  „  nAt  „  nAt  (n+l)At  _(n+l)At 

Ujd*  -  -_S_  rn(j^l)d*  I  n(j-l)d*  n(j+l)d*  ~  n(j-l)d* 
2  1  2d*  2d* 


] 


(n+l)At  nAt  „  nAt  „  nAt  „(n+l)At  „(n+l)At 

njd*  ~  njd«  _  -H  ,u(Jtl)d«  ~  u(]-l)d«  .  “(J+l)d«  ~  u(J-l)d«. 

2  1 9M*  9A* 1 


At 


2d* 


(5.1) 


(5.2) 


Substituting  solutions  (3.3)  and  (3.4)  into  these  equations  yields  a  radian 
frequency 

w*  •  + -jfr  tan*1  SqH  sinfkd*)]  .  (5.3) 


At  l2d* 

The  phase  and  group  speeds  are  then 
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and 


C*  =  ^*  =  +  -j^jr  tan  *  *^H*  sin(kd*)] 

CJ  =  af*  "  ±  ^  cos(kd*)  [1  +  (fj*  ^H)2sin2(kd*)]_1  . 


(5.4) 

(5.5) 


6.  Model  Runs 

We  performed  numerical  experiments  to  test  the  application  of  the  analyses 
in  sections  3  and  4  to  models  with  additional  physical  processes  covering  limited 
domains.  We  used  a  barotropic,  free-surface,  explicit,  primitive  equation  model 
(PL1FBE)  to  compare  with  the  time-centered  explicit  analysis.  We  used  a  baro¬ 
tropic,  free-surface,  time -centered,  implicit,  primitive  equation  model  (PL1FBI) 
to  compare  with  the  time-centered  Implicit  analysis.  Hurlburt  and  Thompson  (1980) 
describe  this  second  model  more  fully. 

Both  models  computed  on  1000  x  1000  x  4  km  basins,  with  d  =  20  km.  For 
each  tlmestep  we  Investigated  the  effect  of  varying  the  wave  resolution  by 
Initializing  the  basin  with  0.5,  1,  2.5,  5  and  10  waves  in  the  basin.  We  ran 
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both  models  for  At.  •  30  sec  and  ran  PL1FBI  additionally  at  At  =  178.6  sec  and 
505.1  sec. 


7.  Discussion 

The  one-dimensional  Courant-Friedrich-Lewy  (C.F.L.)  stability  criterion 
(neglecting  advection)  requires  that  the  gravity  wave  speed  of  the  fluid  remain 
less  than  that  which  would  traverse  the  grid  mesh  distance  d  in  timestep  At, 
that  is  /gi?  £  1.  We  define  the  Courant  number  for  a  one-dimensional  mesh  as 

c0(1)=^5h^  (7.D 

and  for  two  dimensions 

C0(2)=/*iH^  (7.2) 

This  number  gives  the  amount  by  which  we  either  exceed  or  fall  within  the 
theoretical  C.F.L.  criterion.  Table  1  summarizes  the  relevant  phase  and  group 
speed  formula  in  terms  of  CQ(2)  and  kd*. 

Figures  1  and  2  stand  as  checks  on  the  analysis.  Figure  1  compares  the 
analysis  in  section  5  with  values  taken  from  a  plot  resulting  from  a  similar 
analysis  in  two  dimensions  performed  by  Mesinger  and  Arakawa  (1976).  Figure  2 
compares  the  analysis  in  sections  3  and  4  with  results  calculated  from  the  expli¬ 
cit  and  centered  in  time  implicit  models  discussed  in  section  6.  First,  note 
the  agreement  between  the  different  results,  which  indicates  the  usefulness  of 
the  analysis.  Second,  note  the  effects  of  spatial  and  temporal  resolution. 

Figure  3  summarizes  these  effects.  As  expected  for  all  schemes,  gravity 
waves  slow  as  wave  resolution  by  the  grid  spacing,  d,  decreases.  For  the 
implicit  schemes,  gravity  waves  slow  as  Courant  number  Increases.  At  sub-C.F.L. 
conditions,  the  explicit  scheme  reproduces  theoretical  phase  speed  better  than 
the  trapezoidal  implicit  scheme  which,  in  turn,  reproduces  wave  speed  better  than 
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CA/gH 


C/VgH  VS  kd* 


Figure  1.  Ratio  of  numerical  to  theoretical  shallow  water  gravity  wave,  phase 
velocity,  c*/*W  ,  vs.  kd*  from  trapezoidal  implicit  analysis  and  from  plot  of 
similar  two-dimensional  analysis  by  Meslnger  and  Arakawa  (1976).  wavelength 


cV5h 


C*/y/gH  VS  kd* 


Figure  2.  Same  as  Figure  1  for  explicit  and  centered  Implicit  analyses  and  for 
comparable  explicit  and  Implicit  numerical  model  runs. 


Figure  1  for  explicit, 
t  C0(2)  =  0.4,  and  for 


trapezr 
the  im| 


implicit,  and  centered 
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Figure  3.  Same  as 
implicit  analyses  a 
C0(2)  =  7.1. 


»il«  n 


.  m  r« 


the  centered  implicit  method.  At  super-C.F.L.  conditions  the  trapezoidal  scheme 
continues  to  outperform  the  centered  implicit  scheme. 

For  a  quantitative  feel  for  these  comparisons,  consider  a  wave  resolved  by 
10  d  (a  rule-of-thumb  minimum  for  adequate  spatial  resolution).  At  CQ(2)  =  0.42, 
the  explicit  scheme  produces  99%  of  the  phase  speed,  the  trapezoidal  implicit 
scheme  98%,  and  the  centered  implicit  scheme  97.5%.  Exceeding  C.F.L.  results 
in  rapid  deterioration  of  calculated  phase  speed.  At  CQ(2)  =  2.5,  the  trape¬ 
zoidal  implicit  scheme  yields  90%  of  the  theoretical  speed  while  the  centered 
scheme  yields  only  75%.  At  CQ(2)  =  7.07,  the  trapezoidal  method  gives  only  64% 
while  the  centered  implicit  scheme  results  in  only  40%  of  the  actual  value. 

For  the  group  velocity,  the  statement  given  for  the  numerical  effects  on 
phase  velocity  also  apply  qualitatively  to  the  group  velocity.  Quantitatively, 
the  effect  of  decreased  spatial  and  temporal  resolution  is  more  drastic  (Fig.  4). 
For  C0(2)  =  0.42,  and  lOd/wave,  the  explicit,  trapezoidal  implicit,  and 
centered  implicit  schemes  yield  97,  94,  and  92%  of  theoretical  group  speed, 
respectively.  For  C0(2)  =  2.5,  the  trapezoidal  and  centered  implicit  methods 
give  73  and  43%,  respectively.  For  C0(2)  =  7.07,  the  trapezoidal  and  centered 
schemes  yield  28  and  9%,  respectively. 
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8.  Conclusions 

Consider  now  the  task  of  choosing  a  numerical  scheme  to  simulate  a  situation 
where  gravity  waves  may  be  important.  This  situation  might  include  storm  surqe 
and  tsunami  predictions  or  just  forecasts  of  transient,  wind-driven  circulations 
in  ?  coastal  region.  How  should  we  choose  between  the  three  schemes  analyzed? 

If  the  horizontal  resolution  only  minimally  resolves  the  wave  form  and  a  sub-C.F.L. 
timestep  is  affordable,  choose  the  explicit  scheme.  If  we  require  C0(2)  >  1, 
perhaps  for  economic  reasons,  use  the  trapezoidal  implicit  scheme,  but  keep  in 
mind  its  limitations.  If  spatial  resolution  allows  over  50d  per  wave  and  super- 
C.F.L.  timestep  is  required,  either  of  the  implicit  techniques  can  be  used. 

In  summary,  this  report  demonstrates  the  importance  of  and  relative  ease  of 
analysis  in  determining  the  effect  of  time-differencing  schemes  on  shallow  water 
gravity  wave  simulations. 
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decreasing  temporal  resolution.  For  At  JgFVd  <  I,  the  slowing  increases  as  we 
switch  from  explicit,  to  trapezoidal  implicit,  to  centered  implicit  methods.  In 
this  formula,  /gft  =  theoretical  shallow  water  gravity  wave  speed.  At  =  time 
increment,  and  d  =  smallest  space  increment  between  like  variable  grid  points. 
For  At  /gl-j/d  >  1,  the  trapezoidal  implicit  scheme  again  outperforms  the  centered 
implicit  scheme. 
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